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Abstract. The Relativistic Motion Integrator (RMI) consists in integrating numer- 
ically the EXACT relativistic equations of motion, for a given metric (corresponding 
to a gravitational field at first post- Newtonian order or higher) , instead of Newtonian 
equations plus relativistic corrections. 

The aim of the present paper is to validate the method, and to illustrate how RMI 
can be used for space missions to produce relativistic ephemerides of test-bodies (or 
satellites). Indeed, nowadays, relativistic effects have to be taken into account, and 
comparing a RMI model with a classical keplerian one helps to quantify such effects. 

LISA is a relevant example to use RMI. This mission is an interferometer formed 
by three spacecraft which aims at the detection of gravitational waves. A precise orbit 
model for the LISA spacecraft is needed not only for the sake of satellite ephemerides 
but also to compute the photon flight time in laser links between spacecraft, required 
in LISA data pre-processing in order to reach the gravitational wave detection level. 

Relativistic effects in LISA orbit model needed to be considered and quantified. 
Using RMI, we show that the numerical classical model for LISA orbits in the 
gravitational field of a non-rotating spherical Sun without planets can be wrong, 
with respect to the numerical relativisitic version of the same model, by as much as 
about ten kilometers in radial distance during a year and up to about 60 kilometers 
in along track distance after a year... with consequences on estimated photon flight 
times. 

We validated RMI numerical results (using a metric following the International 
Astronomical Union -lAU- 2000 resolutions) with an analytical developpement (up 
to first order in eccentricity and up to first order in GM/(? , where G is Newton's 
constant, M, the solar mass and c the speed of light in vacuum). 

Finally, the RMI relativistic numerical approach is soon more efficient than 
the analytical development. Moreover, RMI extends to other cases (planetocentric, 
instead of barycentric) and can be applied to other space missions. 
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1. Introduction 



Due to considerable increase of the accuracy level in modern space 
missions in the recent years, or expected in close-future missions, rela- 
tivistic gravitational effects must be considered when computing space- 
craft ephemerides. 

Indeed, the Schwarzschild radius {2GM /c?) of the Earth is of the order 
of one centimeter; while that of the Sun is of the order three kilometers. 
The first corresponds to the order of magnitude of the precision in 
current space geodesy; while the second, to the precision requested in 
some future space mission such as in LISA (2015). The rclativistic 
Lense-Thirring effect has already been partially detected with LA- 
GEOS Earth orbiting satellites (Ciufolini & Pavhs 2004). Numerical 
integrations in the post-Newtonian approximation versus Newtonian 
ones have shown the relevance of relativistic effects in the orbit of the 
future GAIA mission (2011) (Khoner 2005), since the GAIA spacecraft 
must be controlled with an accuracy of 0.6mm/s. Owing to the above 
motivations, the present work is dedicated to a numerical relativistic 
model for a generic space mission. 

The method RMI (Relativistic Motion Integrator), a fully consistent 
general relativistic approach (Pireaux et al. 2005), (Pireaux et al. 2006) 
consists in integrating numerically the EXACT relativistic equations 
of motion for a given metric. The advantages of the method are the 
following. All relevant relativistic effects are taken into account if a 
gravitational metric adapted to the precision of measurements is cho- 
sen. The approach is relativistically consistent, and safer than adding 
relativistic corrections by hand to a computation first developed in 
a Newtonian framework. The RMI approach natively contains all the 
gravitational classical and relativistic effects at the corresponding order 
of the metric, including all the couplings between these effects at the 
corresponding order with respect to the metric chosen. This is a serious 
advantage over a Newtonian-plus-relativistic-corrections approach such 
as is implemented in commonly used orbit determination softwares. 
These perturbation approaches become more and more questionable 
as the requested precision increases, requiring a larger number of rela- 
tivistic effects to be taken into account. RMI could help to point out 
deficiencies in common softwares. 

The standard approach to integrate the relativistic differential equa- 
tions of motion are the Einstein-Infeld-Hoffmann (EIH) equations of 
motion (see (Brumberg 1992), (Brumberg 2004), (Brumberg 2007), 
(Damour et al. 1991), (Damour et al. 1992), (Damour et al. 1993), 
(Damour et al. 1994), (Moyer 2000) and references therein). EIH equa- 
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tions are an analytical first order post-Newtonian (IPN) development 
of the exact relativistic equations of motion. The advantage of the 
RMI method over the standard integration of EIH equations is that 
RMI straightaway numerically integrates the equations of motion for a 
chosen metric provided at a given PN order (whether IPN or higher). 
Hence, if, according to new lAU (or else) resolutions, a more appropri- 
ate/precise metric than the present Barycentric Coordinate Reference 
System (BCRS) metric for the motion in the Solar System or Geocen- 
tric Coordinate Reference System (GCRS) metric for planetocentric 
motion (with the Earth as central body) is recommended, RMI can 
straight away use that new metric... without the need to recalculate 
and implement new analytical developments. Only the metric module 
in the RMI software changes. Indeed, separate modules in the RMI 
numerical method also allow easy adaptations and updates for a given 
mission (number of plane containing satellites, number of satellites per 
plane, initial conditions -positions and velocities-), central body pa- 
rameters (mass multipole development of the gravitational potential, 
spin), planetary ephemerides, lAU recommendations (metric, space- 
time transformations)... while keeping the main body of the software 
unchanged. 

When wishing to illustrate how the RMI method can be used in 
space missions, LISA is a good candidate. However, the aim of the 
present paper is not to provide a thorough model of the LISA detector. 
Although some results obtained with RMI for LISA's orbit model are 
relevant for a LISA simulator. 

The LISA (Laser Interferometer Space Antenna) mission (LISA 2000) 
is a space detector of gravitational waves in the [~ 10~^, ~ 10~^] Hz fre- 
quency band. Gravitational waves crossing the LISA quasi equilateral 
triangular constellation are detected through the induced change in the 
station inter-distances. The latter also depend on time, mainly due to 
the gravitational field of the Sun (Chauvineau et al. 2005) around which 
LISA rotates, 20 degrees behind the Earth, and to that corresponding 
to planets; what we call "geometry effects". 

"Noise effects" in LISA are orders of magnitude larger than "gravita- 
tional wave source effects". In order to reach the gravitational wave 
detection level, a Time Delay Interferometry (TDI) method (see (Dhu- 
randhar et al. 2002), (Estabrook et al. 2000)) must be applied to get rid 
of (most of) the laser frequency noise and optical bench noise. The TDI 
method consists in combining numerically data fluxes at the stations 
(rather than combining the laser beams physically) with an appropriate 
delay. Hence, the so-called TDI observables are symmetrized combina- 
tions of the different laser links with appropriate delays (combination 
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of photon-flight time between two stations which correspond to station 

inter-distances) that cancel (almost all) the laser frequency noise and 
optical bench noise. The TDI method is the crucial pre-processing of 
LISA data, before even trying a given strategy to detect any gravita- 
tional wave signal. 

Therefore, in order to validate the new TDI technique and since a 
laboratory replica of the complex LISA mission is not totally achiev- 
able, the performance of LISA TDI can only be studied with com- 
puter simulations of the different processes involved. Such is the aim 
of the LISACodc software [(Pctiteau et al. 2008)] developed by the 
LISAFrance group [(LIS A- France 2009)], or of other simulators in the 
USA [(Vallisneri 2005), (Cornish et al. 2004)]. Among the processes to 
be implemented in a LISA simulator, the orbit model of the spacecraft, 
providing positions, velocities and intcrdistances of spacecraft needed 
for TDI, is the subject of the present paper. 

Relativistic effects in LISA needed to be considered and quantified. 
In the framework of the LISA mission, in articles (Chauvineau et al. 
2005) (see references therein for a generic approach) and (Pireaux 
2007), the photon flight time problem, also sometimes referred to as 
time transfer, and proper time scales of LISA spacecraft are tackled 
using a consistent general relativistic approach. However, the orbit 
model used to compute the initial positions and velocities of LISA 
spacecraft at emission time needed in the time transfer simulation or 
in propcr-versus-coordinate time transformations is classical. 
And so is it still presently the case too in the TDI simulators named 
Synthetic LISA (Vallisneri 2005), LISA Simulator (Cornish et al. 2004) 
and LISACode (Petiteau et al. 2008). 

In the preliminary optimal orbit design for LISA used by Hughes (Hughes 
2005), LISA'S orbit model is also purely classical (in presence of a spher- 
ical non-rotating Sun with planets). The author looks for the optimal 
set of orbital inclinations, eccentricities, semi-major axis, longitude of 
the ascending nodes, arguments of perigee and initial mean anoma- 
lies (afe, Cfc, ifc, fifc, Wfc, Mfco) of LISA spacecraft (A; = 1,2,3) in order to 
minimize LISA's arm flexing according to certain optimization criteria. 

In the present article, we use RMI (assuming no non-gravitational 
forces for LISA spacecraft motion) to quantify the errors implied when 
a classical orbit model is adopted for LISA instead of a general rel- 
ativistic one for the same initial conditions (Barycentric Coordinate 
Reference System -BCRS- positions and velocities of spacecraft). We 
first investigate the case of a classical circular orbit of reference around 
a spherical non-rotating Sun without planets, which we call the circular 
spherical symmetric case. We then extend to eccentric orbits and name 
this case the eccentric spherical symmetric case (more specifically for 
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LISA, e ~ 0.0096). 

Our numerical estimate of relativistic versus Keplerian orbit model 
for LISA with e = 0.0096 shows that the difference between predicted 
barycentric relativistic and classical radial distance reaches up to about 
8 — 9 km during a one-year mission and that the along track difference 
in orbits is about 54 — 59 km after one year (i.e. after one classical pe- 
riod), according to the spacecraft considered, in the eccentric spherical 
symmetric case. 

The relativistic versus classical modelling of LISA's orbit has reper- 
cussions on the flexing of LISA intcrfcrometric arms, the so-called 
breathing of the LISA constellation around its nominal arm-length 
value L = 5 • 10^ m. We show that a relativistic orbit model is relevant 
when studying photon time transfer needed in the TDI method; more 
specifically becaiisc the zcroth order is but the spacecraft inter-distance 
divided by the speed of light. 

Since LISA eccentricity is small and because TDI and classical orbit 
models for LISA used by the Mock LISA Data Challenge (MLDC) 
(Arnaud et al. 2007) task force have been developped using first-order in 
eccentricity approximations, we provide a relativistic analytical check: 
a development up to first order in eccentricity and up to first order 
in GM/c^, where G is Newton's constant, M, the solar mass and c 
the speed of light in vacuum) circular or eccentric spherical symmetric 
cases. 

For the circular spherical symmetric case, the analytical development 
up to first order in e and GM/c^ (equations (11) and (12) ) leads to 
small residuals (about 1 cm in x-y-positions or along track distance and 
a few millimeters in radius) with respect to the RMI numerical rela- 
tivistic model for LISA. However, in the eccentric spherical symmetric 
case, even for a small eccentricity such as LISA's (e ~ 0.0096), the 
corresponding residuals arc non negligible (reaching up to about 85 m 
in along track distance) due to the and higher terms neglected in the 
analytical development; whereas RMI implicitly contains all order in e. 
Hence, the analytical development is soon surpassed by the numerical 
relativistic approach of RMI. This remark is even more relevant to 
space missions with important eccentricities. 

The RMI method was furthermore validated in reference (Hees Sz Pireaux 
(2009)) (for BepiColombo or MarsNext mission) using a full IPN de- 
velopment. 

The present paper is organized as follows. In Section 2, we recall 

the classical orbit model for LISA around a spherical symmetric Sun, 
which is to be our trial example for RMI along this paper. 
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In Section 3, we summarize the RMI relativistic numerical method and 
apply it to LISA with the appropriate initial conditions corresponding 
to the classical model. The numerical results obtained for LISA are 
then discussed. 

In Section 4, we provide an analytical developpement (up to first order 
in eccentricity and in GM/c?) to check RMI. 

Finally, in Section 5, we conclude on the relevance of the RMI approach 
and on the main results obtained for LISA. 

The annex A discusses the numerical accuracy of the RMI method for 
LISA. 

We adopt Einstein's summation convention on repeated indices. 
Latin indices are for space coordinates, such as Z = 1,2, 3; while Greek 
indices arc for space-time coordinates, such as a = 0, 1, 2, 3 with 

^a=0,l,2,3 ^ (^c-t,X,y,z). 



2. LISA classical orbit model in the spherical symmetric 

case 

Presently, within simulators testing LISA TDI (in the framework of 
the LMDC (Arnaud et al. 2007)), the following simplifications relative 
to LISA orbits are assumed. Each spacecraft follows perfectly a free- 
falling test mass that is itself perfectly shielded from non-gravitational 
forces and feels no constraints (for simplicity, one test-mass per space- 
craft is modeled). As the gravitational field is concerned, solely a spher- 
ical non-rotating Sun is considered. The orbit model is classical. 
In present LISA simulators for TDI, departures from the above as- 
sumptions on orbits are presently considered as part of the noise bud- 
get in TDI: among residual laser frequency and optical bench noises, 
scattered-light noise, detector shot noise, laser-beam pointing instabil- 
ity, acceleration noise, inertial noise and others (as specified in Table 1 
of reference (Petiteau et al. 2008)). 

For such a classical orbit model for the three LISA spacecraft k = 
1, 2, 3, in the BCRS, as in (Dhurandhar et al. 2005), the barycentric co- 
ordinates {xk,yk, Zk), for arbitrary initial conditions, can be rewritten in 
terms of rotated Keplerian ellipses {xeii k-. Veil ki Zeii k) with eccentricity 
e ~ 0.0096 as 

f Xk\ / Xell k \ 

[Vk \ = 5i-M Veil k (1) 
\Zk / \ Zell k J 

with 
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/ Xell k \ 




1 Hell k 1 




\ Zell k ) 





a (cos '^k — e) 
aVl — sin '^^ 


where a, e, z and (j are the common semi-major axis, eccentricity, 
inchnation and argument of the periaster of the three spacecraft orbits, 
respectively; and 

= ( Pi P2 P3 ) , 

where the columns of the inverse rotation matrix are given by 

+ cos fifc COS uj — sin Q.k sin u cos i 
pi=\+ sin Q,k cos ui + cos sin uj cos i 
+ sin u} sin i 

— cos Jlfe sin a; — sin $7^ cos a; cos i 
p2 = I — sin rife sin u + cos 0^ cos u) cos i 
+ cos a; sin i 

+ sin Q,k sin i 
p3 = I — cos r^fc sin i 
+ cos i 

Indeed, we start from a slightly different hypothesis with respect to 
Hughes' (Hughes 2005). We take common {a,e,i) for the three space- 
craft with optimal e, i in order to minimize LISA's arm flexing in 
agreement with reference (Nayak et al. 2006): 
a = 1 AJJ^ 

, / ^ sin I 

where u = f + 1^ is the optimal inclination of the LISA triangle on the 
ecliptic and L = 5-10^ m is the average interferometric arm- length. The 
longitude of the ascending node, $7^, is particular to a given spacecraft 
k and is given in terms of that of the first one with a phase shift i?^: 

nk = ni-'dk withi?fe = -2(/c-i)|. 

The time parametrization of the orbits is given by the equation of the 
eccentric anomaly '^k of each spacecraft, 

*fe-esin*fe = Mfc , (2) 

with the mean anomaly 

Mfe = ^ (t - to) + Mko 
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in terms of the orbital period, T, and the mean anomaly of spacecraft 
k at initial time to, that is M^g = Mk{t = to)- 

Mean anomalies are related to that of the first spacecraft through the 
phase shift: 

Mfc = Ml + . 

BCRS position and eccentric anomaly equations used in (Chau- 
vineau et al. 2005) correspond to particular initial conditions {to = 
0, Lo = 3-7r/2, 0,1 = 37r/2, Miq = 0) without any planets (which means 
that both the initial time, to, and the initial mean anomaly of the first 
spacecraft, Mio, are completely arbitrary in that case). 
We also recall that the time when spacecraft k is at perihelion is given 

by 

Mko 



with the mean motion n = 27r/r = ^GMj from Kepler's 3rd law. 



3. Numerical native relativistic orbit model 

3.1. Exact relativistic equations of motion 

In General Relativity, the motion of a spacecraft is described by the 
relativistic equation of motion. 



1^ 



-pa 



dx^ dx"' 
dr dr 



dx" dx^ 
dr dr 



(3) 



where Kp is a quadri- "force" encoding non-gravitational forces; r, the 
proper time aboard the considered spacecraft; and T^^, Christoffel 
symbols with respect to the metric. The relation between covariant 
and contravariant metric components being 



9 ■ 9M = K- 



(4) 



The four equations in (3) are redundant because of the normalization 
of the quadrivelocity. 

In the case of LISA, assuming only one shielded test-mass per satel- 
lite, each satellite follows a geodesic motion, that is Kp = 0. Combining 
equations in (3) , we can remove the proper time variable to rewrite the 
set of relativistic equations as 



dh 



df^ 



dx^ dx^ 
~dt dt 



(5) 
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3.2. Relativistic Motion Integrator (RMI) method applied 
TO LISA 

The Relativistic Motion Integrator (RMI) method (Pireaux et al. 2005), 
(Pircaux ct al. 2006), consists in integrating numerically the exact rel- 
ativistic equations of motion (3) for a given metric. 
The numerical accuracy and stability of the RMI method for the LISA 
mission is validated in Annex A. 

When using the RMI method for LISA, rotating around the Sun, 
the appropriate metric is the BCRS metric recommended by the lAU, 
International Astronomical Union, 2000 resolutions (see (Soffel et al. 
2003) and references therein) and the corresponding isotropic coordi- 
nates. The BCRS lAU 2000 metric neglects only terms at order 
and above in g^^ or g^'^; and at order and above in (7'™. The lAU 
2000 resolutions have been adopted in 2000 so to take into account the 
best precision of present and next future space experiments. That is 
experiments involving (or which can be translated in terms of) clocks, 
with accuracies better than a few parts in 10-*^^ in fractional frequency 
and stabitities better than about CTyi^r) = 1-10~^^t~^/^ (Allan standard 
deviation), located at distances as close as 0.25 A.U. from the Sun 
(Soffel et al. 2003). 

Note that most NASA and ESA space missions are modeled accord- 
ing to the EIH equations and corresponding relativistic algorithms 
described e.g. in (Moyer 2000). Unfortunately, reference (Moyer 2000) 
was published around October 2000 and thus does not take into account 
the latest IAU2000 resolutions, published later. (Moyer 2000) refers to 
lERS 1997 resolutions at the latest. 

3.3. LISA INITIAL CONDITIONS 

We shall use the subscript *ci for classical quantities and *rei for 
the relativistic ones. 

In our problem of comparing relativistic and classical LISA ephemerides 
(E), we chose to take the same initial conditions (IC) in terms of 
coordinate positions and velocities of spacecraft A; = 1, 2, 3 for both the 
relativistic and the classical orbits. Indeed, we could have chosen to 
speak in terms of same energy and momentum, but this does not reflect 
the way the actual space mission will be planned and this does not 
easily provide insight in terms of what is the error in predicted position 
and velocities. Hence, initial conditions of the relativistic model will 
be those BCRS {xk,yk, z^; dxk/ dt, dy^/ dt, dz^/ dt) obtained by setting 
t = to in the classical equations (1) and (2). 

Note that this choice is not restrictive since, if the classical and rela- 
tivistic IC differ. 
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Erel{ICrel) — Ecl{ICcl) — [Erel{ICrel) — Ecl{ICrel)] 

+ [E,i{ICrel) - E,i{ICd)] 

and the second r.h.s. term in the above equation, not discussed in this 
paper, is but a classical problem. 

The eccentricity of the numerical ephemerides for the eccentric case is 
that corresponding to LISA spacecraft, e ~ 0.0096. In our numerical 
simulation, we arbitrarily further chose to = 0, u = 37r/2, J7i = 37r/2 
and Mio = in agreement with the initial conditions of paper (Chau- 
vineau et al. 2005). 

Let us point out that this analysis could have been applied to Hughes' 
initial conditions (Hughes 2005). 

3.4. Discussing numerical results for LISA in the spherical 

SYMMETRIC CASE 

The spherical symmetric model for LISA corresponds to a classical 

orbit of reference around a spherical non-rotating Sun without planets. 
Owing to this symmetry, the value of the inclination i is irrelevant in 
order to compare relativistic versus classical ephemerides generated for 
LISA. Hence, wc used the classical method without planets, described 
in Section 2, with i set to and e set to either (circular case) or 
0.0096 (eccentric case), to produce a numerical classical ephemeris 
for LISA {xk,yk,Zk,dxk/dt,dyk/dt,dzk/dt)t. We then used the RMI 
method, described in the above Paragraphs 3.2 and 3.3 with identical 
initial conditions, to produce a corresponding relativistic numerical 
ephemeris. We then used those two ephemerides, recorded as a function 
of BCRS time, to plot (relativistic - classical) quantities as a function 
of BCRS time every day during 365 days (~ T = 27r/ra) such as in 
Section 6. 

3.4.1. Circular classical reference orbit case: 

From Figures 1 and 2, we found that the difference between predicted 
barycentric relativistic and classical x-y-positions reaches up to a max- 
imum of about 51 — 56 km during a one-year mission. 
When speaking in terms of a difference in radial or along track distance 
between numerical relativistic and classical orbits, the above cited re- 
sults translate into Figures 3 and 4, respectively. We computed that 
the maximum difference in radius is about 8.9 km while the along 
track difference in orbits after one classical period is about 56 km for 
this circular spherical symmetric case. 

The spacecraft is ahead on the classical orbit with respect to the rela- 
tivistic one. 

We see from Figure 3 that, having adopted a circular classical orbit of 
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reference, the corresponding relativistic orbit is non-circular. 
The difference in velocity components along the x- or y-BCRS axis as 
a function of time obtained are given in Figures 5 and 6, respectively. 
The difference between predicted barycentric relativistic and classical 
x-y-velocities reaches up to a maximum of about 0.007 — 0.010 m/s 
during a one- year mission. This agrees with the order of magnitude for 
the difference in position over one year. 

3.4.2. Eccentric (e = 0.0096j classical reference orbit case: 
Prom Figure 7, we see that the maximum difference in radius between 
numerical relativistic and classical orbits is about 8 — 9 km, according 
to the spacecraft considered. From Figure 8, we see that the along track 
difference in orbits after one classical period reaches about —59 or —54 
km, according to the spacecraft considered, for this eccentric spherical 
symmetric case. 



3.4.3. LISA 's arm flexing and photon time transfer: 
Assuming e = 0.0096 and using the numeric relativistic cphemerides 
for LISA spacecraft obtained with the RMI method or that obtained 
with a classical method, we can compute the interferometric-arm length 
Ljk, that is the interdistance between spacecraft j and k. Over a year, 
LISA constellation shows some breathing or triangle flexing: the relative 
position of spacecraft varies as a function of time. It is interesting 
to see that, for the uninclined {i = 0) eccentric spherical symmetric 
model, the classical approach is wrong by as much as about 4 km 
over a one-year mission. However, the true mission has an inclination 
i such as to minimize the breathing (Nayak et al. 2006). Figure 15 
illustrates LISA breathing in the inclined (with the appropriate i given 
in Section 2) eccentric spherical symmetric case. In that realistic model, 
the classical approach is wrong by as much as about 3 km over year 
of mission, as shown by the residuals (relativistic - classical) relative 
positions of spacecraft in Figure 16. This error translates into a missing 
~ 1 • 10~^ s at zeroth order in GM/{a (?) oc v^jc? in photon time 
(0) 

transfer (tjk= Lji^/c) after a year. We recall that, in paper (Chauvineau 
et al. 2005) where the time transfer of photons between LISA spacecraft 
was studied for a classical LISA orbit, the zeroth order amounted to 
about 16.7 s (5 • 10® km/c, that is the nominal interferometric arm- 
length, L, traveled at the speed of light) with a flexing amplitude of 
about 0.16 s (48000 km/c); the half order amounted to about 3 • 10~^ 
s (960 km/c); and the first order was less than about 1 • 10~^ s (< 30 
m/c). Hence, we understand the relevance of relativistic orbit model in 
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the TDI approach, for a coherent modeUing of the mission over a few 
months. 



4. An analytical development in eccentricity to check the 
numerical relativistic versus classical orbit model 



Let us find an analytical check of the (relativistic - classical) numerical 
integration in the eccentric spherical symmetric case, up to first order 
in e and GM/c^. At the post-Newtonian level the solution is known in 
terms of osculating elements or other representations (e.g. in (Brumberg 
1991) or Annex 2 in (Soffel 1989)), valid for any eccentricity. However, 
those are implicit solutions (for the radial distance and polar angle) and 
a further development in eccentricity would be relevant to the LISA 
mission. Indeed, in present LISA literature, orbits and Time Delay In- 
terferometry (TDI) are considered at different levels of approximation, 
based on a (classical) development in terms of the small eccentricity of 
the LISA mission (an orbit development at a first-order in eccentricity 
is further assumed by (Arnaud et al. 2007)). For example, to be ideally 
a 100 percent efficient in removing laser frequency noise and optical 
bench noises, the TDI combinations from 1st generation TDI algebra 
assume symmetric and constant (in time) photon propagation time 
between two LISA spacecraft. This is met only by a rigid motionless 
constellation model. Hence the need for a 1.5th TDI generation algebra, 
this time relaxing the symmetry on time-delays. The latter TDI as- 
sumptions being met by modeling the constellation as rotating around 
its center of mass, and around the Sun (without any planet present) 
in a Keplerian motion at first order in eccentricity. Deviations from 
this 1st order in eccentricity Keplerian model lead to residual laser 
frequency and optical bench noise in the TDI combinations, which 
need to be quantified. Consequently, the explicit general relativistic 
solution provided in this section as a development at IPN and first 
order in eccentricity is useful for the sake of comparison with existing 
LISA classical models. Our analytical development provides the explicit 
{Srk = rk rel - rk cl, S9k = Ok rel - Gk cl) relativistic upgrade to the 
Keplerian 1st order in eccentricity orbit model for LISA such as used 
by the LMDC (Arnaud et al. 2007). 

To proceed, we first develop the geodesic equation of motion (5) 
up to the corresponding order in GM/(? in the BCRS. Writing = 
4ei we find 
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where = dx^ /dt is the velocity of spacecraft at time t in the BCRS. 
Using the analytical developments of ChristofFel symbols in the BCRS 
at the corresponding order, we can write the difference between the 
relativistic and classical orbit accelerations d^e^ /dt^ as 
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(6) 



where r is the coordinate radial distance relative to the Sun at time t 
in the BCRS and (s) means that the term considered is of order s in 
GM/c^. 

Since we consider a symmetric gravitational field and are interested in 
the difference between relativistic and classical ephemerides for a given 
satellite, the inclination i is irrelevant. Hence, we choose to work with 
i = 0. The inclined analytical solution can be obtained by a simple 
rotation of the uninclined analytical solution (11, 12). Then of course, 
Zci = Zrei = = 0, as well as the corresponding time derivatives. 
Let us further use the set of polar coordinates (r, 9) with x = rcosO, 
y = rsinO and 2; = to reflect the symmetry of the problem. The 
above set of equations (6) becomes 
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(8) 

where 69 = 9^,.^ — 9(.i, 6r = rrei — Td and *= d* jdt. 

Using Kepler's orbital motion equations (r^; = a(l — e^)/(l + ecos9d), 

r^l 9= ^GM a (1 - e^), v^i = ^JGM {2/rd - 1/a)), we can check that 
equations (7) and (8) lead to two first integrals of the motion: 

rl-5 9+2 Td- 9d -Sr + ^^JgM a (1 - e^) = Q (9) 



Td -Sr + 



•2 GM 

^cl +— 2- 



Sr + 9d-S9^ 



GM 



rti rd 



= Dk 
(10) 



Those can be traced back to the relativistic angular momentum and 
energy integral resulting from the spherical symmetry. 
Owing to our choice of identical positions and velocities of spacecraft 
at initial time for both the classical and the relativistic orbit models, 

{69, 6r, 6 9,6 r)t^ = (0, 0, 0, 0). Hence the integration constants are 



Ck 



Dk = 



GM ^GM a (1 -e^) 



GM 



fk cl 



3 ^fc cl 
rk cl 



GM 

' k cl 



with Tk cl 0= Tk clito) and Vk d o = Vk dito) given by Kepler's orbital 
equation of motion at initial time with respect to the initial conditions 
of a given spacecraft /c = 1,2 or 3. Equations (9, 10) provide a first 
check of the numerical results of Sections 3.4.1 and 3.4.2 in the spherical 
symmetric approximation. We note that for a circular orbit of reference 

{ud =9d, fd = fd= 0), Ck and the third term of the left-hand-side 
of (9) cancel; while Dk and the fourth term of the left-hand-side of (10) 

cancel... leading to the same identical first integral: 6 1= —2 Ud 6r, 
where 61 = r^ei ■ 59. 

We now develop the differential system ((7) and (8); or, which is 
easier, (9) and (10)), up to first order in e using Kepler's equations of 
motion at first order in e: 
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with *= d* /d{ncit). To find solutions to the above differential system, 
we use the theory of perturbation around nuh eccentricity. We find 
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[1] -{3 + llcos2(nc/tfep)}cos(nc/(t-tfcp)) 

^rjfc = +e—^{ +11 cos {ridtkp) sin {ridtkp) sin {nd{t- tkp)) 
-6 cos {ridtkp) cos^ (nci(i - t^p)) 
+6 sin {ridtkp) sin (nci(t - tkp)) cos (nci(t - tkp)) _ 

(12) 

where [s] means that the term considered is of order s in e. At zeroth 
order in e, those results correspond to the circular classical orbit of 

reference case. 

Expressions (11) and (12) can be easily transposed in terms of (relativis- 
tic - classical) positions {Sxk, Syk) and related (relativistic - classical) 

coordinate velocities {S Xk, 6 Uk) using 



Xk = Vk COS 6k 
Vk = rk sin 9k 



Sxk = cos Ok ■ Svk ~ Vk sin 0k ■ 50k 
Syk = sin 0k ■ Srk + cos 9k ■ 50k 
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4.1. Circular classical reference orbit case: 

Expressions (11) and (12) with e = match perfectly the numerical 
results for the circular spherical symmetric case presented in Section 
3.4.1, up to first order in GMjc?. Residuals between RMI approach 
and this analytical check for the circular spherical symmetric case 
reach about 1 cm in x-y-positions or along track distance and a few 
millimeters in radius (Figures 9 and 10). 

A dimensional analysis leads to an order of magnitude for the difference 
between classical and relativistic barycentric positions of spacecraft of 
about GMj (ac^) -27^ ~10 km for a one year simulation. Our numerical 
native relativistic approach shows that classical modelling can be wrong 
by as much as about 50 km, in terms of barycentric coordinates (x,y,z) 
and along track distance, over one year. It is interesting to point out 
that this is nearly one order of magnitude larger than estimated with a 
dimensional analysis. The numerical results are confirmed by the more 
cautious analytical developpements presented above. 

4.2. Eccentric classical reference orbit case 

Expressions (11) and (12) with orbital elements corresponding to LISA's 
(e = 0.0096) hnt i = match the numerical results for the eccentric 
spherical symmetric case presented in Section 3.4.2, up to first order in 
e and in GM/c^. 

Residuals between the RMI approach and this analytical check at zeroth 
order in e, for the eccentric spherical symmetric case, reach up to 
about +840, ±540 or —800 m in radial distance and about —3600, 
+2400 or +1600 m in along track distance, for spacecraft A; = 1, 2 or 3 
respectively, over a year (Figures 11 and 12). 

When the analytical check for the eccentric spherical symmetric case 
is considered up to first order in e, the residuals reach up to about 
+24, —15 or +14 m in radius and about —85, —25 or +32 m in along 
track distance, for spacecraft = 1,2 or 3 respectively, over a year 
(Figures 13 and 14). Residuals between the RMI numerical analysis 
(implicitly containing all orders in e) and the analytical development 
(up to first order in e, equations (11) and (12) ) are bound to be 
larger for space missions with larger eccentricities than that of LISA's 
{^LISA — 0.0096). This shows the limits of the analytical development 
even for an eccentric model with a simple spherical symmetric gravi- 
tational field. And going to higher orders in e increases the number of 
terms in expressions (11) and (12) drastically, as illustrated by the 0-th 
and 1-st order contributions. Non symmetric cases such as in presence 
of planets, with a central body which is non spherical or has a spin, are 
even much more complex to handle analytically. On the opposite, the 
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RMI approach, which is exact in terms of e, impUcit in terms of spin, 
flattening or planets, via the metric, is very flexible. Indeed, RMI also 
runs when the spherical symmetry is broken, since solar spin, multipolar 
development of the solar mass, point-like planets can be introduced in 
the metric and hence be coherently taken into account in numerical 
ephemerides produced for LISA via that approach. 



5. Conclusions 

The aim of the present paper was to illustrate how the Rela- 
tivistic Motion Integrator (RMI) can be used to provide a relativistic 
numerical satellite (or test-body) propagator for space missions; and to 
quantify relativistic effects when a comparison is made with a classical 
corresponding model. 

As an illustration of RMI and to validate the method, we chose 
the space interferometer LISA, modelled in the Bary centric Coordi- 
nate Reference System (BCRS) in the gravitational field of a spherical 
non-rotating Sun, without planets (the spherical symmetric case). We 
compared the numerical relativistic ephemeris (propagated daily po- 
sitions and velocities of each spacecraft) obtained with RMI to the 
ones obtained with a classical numerical model with identical initial 
conditions in terms of positions and velocities. The (relativistic - clas- 
sical) BCRS position obtained seemed a priori large, up to a few tenth 
kilometers, i.e. more or less 5 or 6 times the estimate obtained from a 
rapid dimensional analysis. 

However, we made a careful analytical analysis: analytical expressions 
(up to first order in GM/c^, with G, Newton's constant, M, the Sun's 
mass and c, the speed of light in vacuum) of two first integrals of 
the problem and an analytical development of (relativistic - classical) 
BCRS along track and radial distances up to first order in eccentricity 
e and in GMjf?. The analytical developments with orbital elements 
corresponding to LISA's confirmed the numerical results obtained and 
validates the RMI approach. The difference between the RMI numer- 
ical approach, based on the exact relativistic equation of motion with 
respect to the BCRS metric (which is up to second order in GM /(? 
in the IAU2000 resolutions) for a spherical non-rotating Sun, and the 
analytical development are of order ■ GM/(?. 

Hence, for LISA, we have shown that, when the classical orbit of refer- 
ence is eccentric with clisa — 0.0096, the difference between relativisti- 
cally and classically modelled radial distance reaches up to a maximum 
of about 8 — 9 km during a one-year mission. After one year (i.e. one 
classical period) , the difference in orbits in terms of radial distance can 
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be as much as about 680 m and along track difference is about 54 — 59 

km according to tfie spacecraft considered. 

Errors in LISA satellite orbit may have consequences when modelling 
LISA'S arm flexing for the sake of interferometry. We showed that a 
relativistic orbit model is relevant when studying photon time transfer 
needed in the TDI method. Using a classical orbit model contributes 
to an error of about 10~^ s (~ 3 km/c) in photon time transfer over 
a year. The TDI method is the crucial pre-processing of LISA data, 
before even trying a given strategy to detect any gravitational wave 
signal. 

Since the orders of magnitude of (a^, e^, ik)k=i,2,3 used in Hughes orbit 
model for LISA's three spacecraft (Hughes 2005) are the same as the 
ones chosen here, the same conclusions will apply in Hughes'case. 
Note that, in the present paper, we did not aim at a complete model of 
the LISA detector, but some of the above result might be interesting 
when building a LISA simulator. 

Our present study also shows that while the analytical development 
soon reaches its limits, the strength of the RMI approach is that it 
also runs when the spherical symmetry is broken [i ^ 0, non-spherical 
Sun, rotating Sun, with planets), cases much more complex to model 
analytically. Indeed, a solar spin or multipolar development of the solar 
mass (solar J2) or point-like planets can be introduced in the metric 
and hence be coherently taken into account in numerical ephemerides 
produced for LISA via the RMI approach. The point is to use a metric 
with a sufficiently high order of development in 1/c^, so as to include all 
the classical and relativistic effects relevant to the precision of the space 
mission considered. The lAU 2000 BCRS metric models coherently, for 
LISA and other space missions, the action of the Sun and planets at a 
relativistic level. 

Finally, the RMI approach can be applied to other space missions, 
whether barycentric or planetocentric. 
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Figure 1. DifTerence between numerical relativistic and classical position 
ephemerides for the LISA mission in the circular spherical symmetric case: x 
barycentric coordinate (5x). 
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Figure 2. Difference between numerical relativistic and classical position 
ephemerides for the LISA mission in the circular spherical symmetric case: y 
barycentric coordinate {5y). 
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Figure 3. DifTerence between numerical relativistic and classical position 
ephemerides for the LISA mission in the circular spherical symmetric case: radial 
barycentric distance (Jr). 
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Figure 4- DifTerence between numerical relativistic and classical position 
ephemerides for the LISA mission in the circular spherical symmetric case: along 
track distance {SI = rrei ■ 59 ^ a ■ 59). 
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X differential velocities (m/s) 




Figure 5. Difference between numerical relativistic and classical velocity 
ephemerides for the LISA mission in the circular spherical symmetric case: velocity 
component along the x barycentric coordinate axis {5 x). 
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Figure 6. Difference between numerical relativistic and classical velocity 
ephemerides for the LISA mission in the circular spherical symmetric case: velocity 

component along the y barycentric coordinate axis {5 V). 
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Figure 7. DifTerence between numerical relativistic and classical position 
ephemerides for the LISA mission in the eccentric {clisa — 0.0096) spherical 
symmetric case: radial barycentric distance {5r). 
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Figure 8. Difference between numerical relativistic and classical position 
ephemerides for the LISA mission in the eccentric {clisa — 0.0096) spherical 
symmetric case: along track distance {51). 
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Figure 9. Residuals between the numerical (relativistic - classical) position 
ephemerides and the corresponding analytical development for the LISA mission 
in the circular (e = 0) spherical symmetric case: radial distance (5r). 
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Figure 10. Residuals between the numerical (relativistic - classical) position 
ephemerides and the corresponding analytical development for the LISA mission 
in the circular (e = 0) spherical symmetric case: along track distance {51). 



PireauxS_REL_vs_CLASSIC_DRBIT0_250609.tex; 6/08/2009; 8:41; p. 23 



24 



S. Pireaux and B. Chauvineau and A. Hees 



residual radius (m) - 0th order in e. - 




spacecraft 1 

spacecraft 2 
— — — spacecraft 3 

Figure 11. Residuals between the numerical (relativistic - classical) position 
ephemerides and the corresponding analytical development at 0th order in e for 
the LISA mission in the eccentric [clisa — 0.0096) spherical symmetric case: radial 
distance {5r). 
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Figure 12. Residuals between the numerical (relativistic - classical) position 
ephemerides and the corresponding analytical development at 0th order in e for 
the LISA mission in the eccentnc {clisa — 0.0096) spherical symmetric case: along 
track distance (SI). 
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Figure 13. Residuals between the numerical (relativistic - classical) position 
ephemerides and the corresponding analytical development up to 1st order in e 
for the LISA mission in the eccentric {e^jsA — 0.0096) spherical symmetric case: 
radial distance {Sr). 



sidual along track distance (m) . oii + isi ordet in c 




SF^acectafl 1 
sfjacectafl 2 
sfjacectafl 3 



Figure 14- Residuals between the numerical (relativistic - classical) position 
ephemerides and the corresponding analytical development up to 1st order in e 
for the LISA mission in the eccentric {clisa — 0.0096) spherical symmetric case: 
along track distance {51). 
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Figure 15. Numerical relativistic modelling of LISA breathing in the eccentric spher- 
ical symmetric case {clisa — 0.0096): relative positions between spacecraft, with 
Ljk the interdistance between spacekraft j, fc = 1, 2, 3 where j 7^ k. 
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Figure 16. Difference between numerical relativistic and classical modelling of LISA 
breathing in the eccentric spherical symmetric case {e^isA — 0.0096): difference in 
relative positions between spacecraft, with Ljk the interdistance between spacekraft 
j, fc = 1, 2, 3 where j 7^ fc. 
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Appendix 

A. Numerical estimate of the ChristofFel Symbols 

Within a numerical integration of the relativistic equations of motion, 
one has to carefully check the numerical accuracy. In this section, we 
show that the numerical errors are smaller than the order of magnitude 
of the relativistic effects. 

In order to integrate equation (3) , we need to evaluate numerically the 
Christoffel Symbols 

r;^v = ^^""^ {9fiu,^, + 9nP,u - 9t,v,p) (13) 



where f^x = and the matrix gf"^ is the inverse of the matrix g^p 
owing to expression (4). We need to evaluate numerically the derivative, 
9nv,(i-, of the metric components. The derivative is computed using an 
estimation of order 4 (Kincaid & Cheney 2002) 



f{x - 2h) - 8/(x -h)+ 8f{x + h)- fix + 2h) 



(14) 



with fix) = Dh{x) + 0{h^) (15) 

As can be seen in Figure 17, one needs to choose the discretisation step 
size, /i, very carefully. For large /i, the discretisation error is important 
(oc /i^) but for small h, the roundoff error increases (oc l/h). 

In order to increase the precision of the derivative, it is usefull to 
derive h^j^i, = g^i, — rj^j,^, where rj^i, is the Minkowsky metric, instead of 
g^i, as it is more stable from a numerical point of view (Figure 17). 

It is also interesting to use Richardson extrapolation (Richardson & 
Gaunt 1927). This requires two estimations of order 4 (-Do,o = Dh{x) 
and -Di^o = -f^h/fe(a^) where A; is a real factor) to construct a new 
estimation of order 8: 

In practice, the factor k is chooscn as 1.5 or 2 and this procedure can 
be iterated starting from Di^q = Df^^f^i to construct the new estimation 

After n steps, Z?n,n is of the order of 0{h^^'^~^^^). Figure 17 ilhistrates 
in the case of LISA how a relative error of order of 10~^^ on the 
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h(m) 

Figure 17. Representation of the relative precision of gtt,x for one point of the LISA 
orbit. The relative precision of gtt,x and htt,x are represented as function of the 
discretisation step, h. The Richardson extrapolation is also represented for a factor 
k = 1.5 (12 iterations are represented). 



derivative of h^y can be reached (in double precision) using Richardson 
extrapolation. This method does not require to start with a very fine 
tuned initial step size h and it is possible to stop the iterations when 
the convergence is sufficient. 
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